% Example on Multiplication Quantization effects on the 
%  First-Order IIR filter:
%   y(n) = Q[a*y(n-1)] + x(n); 
%   x(n) scaled to prevent overflow

close all; clc;

% Example Parameters
B = 12;               % # of fractional bits
N = 100000;            % # of samples
xn = (2*rand(1,N)-1); % Input sequence - Uniform Distribution
a = 0.9;              % Filter parameter
Xm = 1-abs(a);        % Scaling factor

% Local variables
  bM = 7; DbM = 2^bM;    % bin parameter
  BB = 2^B;              % useful factor in quantization
   M = round(DbM/2);     % Half number of bins
bins = [-M+0.5:1:M-0.5]; % Bin values from -M to M
   Q = bins/DbM;         % Normalized bins
 YTN = 2^(-bM);          % Ytick marks interval
 YLM = 4*YTN;            % Yaxis limit

% Quantize the input and the filter coefficients
  xn = round(Xm*xn*BB)/BB; % Scaled Input quantized to B bits
   a = round(a*BB)/BB;     % a quantized to B bits

% Filter output without multiplication quantization
  yn = filter(1,[1,-a],xn);  % output using filter routine

% Filter output with multiplication quantization
   yq = zeros(1,N);               % Initialize quantized output array
yq(1) = xn(1);                    % Calculation of the first sample yq(1)
for I = 2:N;
    A1Y = round(a*yq(I-1)*BB)/BB; % Quantization of a*y(n-1)
  yq(I) = A1Y+xn(I);              % Calculation of the I-th sample yq(I)
end

% Output Error Analysis
   en = yn-yq;                          % Outout error sequence
varyn = var(yn); varen = var(en);       % Signal and noise power   
eemax = max(en); eemin = min(en);       % Maximun and minimum of the error
enmax = max(abs([eemax,eemin]));        % Absolute maximum range of the error
enavg = mean(en); enstd = std(en);      % Mean and std dev of the error
   en = round(en*(2^bM)/(2*enmax)+0.5); % Normalized en (interger between -M & M)
   en = sort([en,-M:1:(M+1)]);          % 
    H = diff(find(diff(en)))-1;         % Error histogram
    H = H/N;                            % Normalized histogram
 Hmax = max(H); Hmin = min(H);          % Max and Min of the normalized histogram

% Output SNRs
SNR_comp = 10*log10(varyn/varen);
SNR_theory = 6.02 + 6.02*B + 20*log10(Xm);

% Plot of the error histogram
Hf_1 = figure('paperUnits','inches',...
    'paperposition',[0,0,6,4],'color',[0,0,0]);
set(Hf_1,'NumberTitle','off','Name','Example 9.6.7');

bar(Q,H); %axis([-0.5,0.5,0,YLM]);
%title(TITLE,'FontSize',16,'FontName','Times', 'Fontweight','Bold');
set(gca,'TickDir','in','FontWeight','Bold'); 
xlabel('Normalized error'); ylabel('Distribution of error');
set(gca,'XTickMode','manual','XTick',[-0.5:0.1:0.5]);
set(gca,'YTickmode','manual','YTick',[0:YTN:YLM]);
LABEL1 = [' SAMPLE SIZE N = ', num2str(N)];
LABEL2 = ['   PARAMETER a = ', num2str(a)];
LABEL3 = ['  SIGNAL POWER = ', num2str(varyn)];
LABEL4 = ['  SNR (THEORY) = ', num2str(SNR_theory)];
LABEL5 = [' ROUNDED T0 B  = ', num2str(B), ' BITS'];
LABEL6 = ['    ERROR MEAN = ', num2str(enavg)];
LABEL7 = ['   NOISE POWER = ', num2str(varen)];
LABEL8 = ['SNR (COMPUTED) = ', num2str(SNR_comp)];
text(-0.47,3.80*YTN,LABEL1,'Fontname','Courier','FontWeight','Bold','Fontsize',10);
text(-0.47,3.65*YTN,LABEL2,'Fontname','Courier','FontWeight','Bold','Fontsize',10);
text(-0.47,3.50*YTN,LABEL3,'Fontname','Courier','FontWeight','Bold','Fontsize',10);
text(-0.47,3.35*YTN,LABEL4,'Fontname','Courier','FontWeight','Bold','Fontsize',10);
text(0.05,3.80*YTN,LABEL5,'Fontname','Courier','FontWeight','Bold','Fontsize',10);
text(0.05,3.65*YTN,LABEL6,'Fontname','Courier','FontWeight','Bold','Fontsize',10);
text(0.05,3.50*YTN,LABEL7,'Fontname','Courier','FontWeight','Bold','Fontsize',10);
text(0.05,3.35*YTN,LABEL8,'Fontname','Courier','FontWeight','Bold','Fontsize',10);